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1.0 INTRODUCTION 


The following is the final report on the development of two 
computer programs r VSAERO-TS and VSAERO-H, for calculating the 
unsteady subsonic aerodynamic characteristics of arbitrarily- 
shaped wings oscillating in pitch. The work is part of an inves- 
tigation of the combined effects of unsteady motion and planform 
shape on the aerodynamic loading of helicopter blades. Four 
blade-tip planforms were evaluated; these include a rectangular, 
a tapered, a swept and an ogee shape. Figure 1. Each blade is 
considered as a semispan wing oscillating in pitch; i.e., blade 
rotation is excluded at this time. 

Experimental data have been measured on the four blades by 
the Institut fur Aeroelastik in the DFVLR 3 x 3 m tunnel in 
Gottingen, West Germany. Correlations between calculated and 
measured unsteady chordwise pressure distributions were presented 
in Reference 1. The present report describes the theoretical 
background to the methods and includes results from parametric 
studies. 

During the correlation study the time-stepping program had 
difficulties with the tapered and ogee planforms and also with 
the tip vortex effects. Modifications have now been installed in 
the VSAERO-TS program to overcome these difficulties and are 
discussed in this report. 

Both programs, VSAERO-TS and VSAERO-H, have a common basis 
(Figure 2) in program VSAERO, which is under continued develop- 
ment (2) through (6) for the analysis of steady, non-linear 
aerodynamic characteristics of arbitrary configurations with 
extensive vortex separations. Program VSAERO-TS is a time-step- 
ping analysis capable of treating large amplitude motions while 
program VSAERO-H uses a harmonic wake and small amplitude as- 
sumptions. VSAERO-H, therefore, is basically similar to 
Geissler's method (7); the main difference is that VSAERO-H uses 
an internal Dirichlet boundary condition to establish a unique 
relationship between the source and doublet singularity distribu- 
tions. 


The basis of these computer programs is a surface singular- 
ity panel method which uses quadrilateral panels on which doublet 
and source singularities are distributed in a piecewise constant 
form. The panel source values are directly determined by the 
external Neumann boundary conditions controlling the normal com- 
ponent of the local resultant flow; the doublet values are solved 
after imposing the internal Dirichlet boundary condition of zero 
perturbation potential at the center (underside) of all the 
panels simultaneously. Surface perturbation velocities are ob- 
tained from the gradient of the doublet solution, while field 
velocities are obtained by direction summation of all singularity 
contributions. The details of the mathematical formulation are 
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Figure 1. Planforms for Unsteady Blade-Tip Study. 
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Figure 2. VSAERO Unsteady Program Development. 
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presented in Section 3, while the numerical procedure for ob- 
taining the unsteady pressure distribution, forces and moments is 
presented in Section 4. 

In Section 5, the viscous/potential flow iterative scheme is 
discussed. This capability comes from the basic VSAERO program 
(Figure 2) and is applied here in VSAERO-TS at the mean angle of 
attack for the purpose of establishing the extent of tip-edge 
separation. A typical calculation for steady flow predicting the 
tip flow is presented for the rectangular tip. 

In Section 6, a discussion on the convergence characteris- 
tics of the VSAERO-H wake integrals is presented. The two doub- 
let distribution models for representing the wake are discussed. 

In Section 7, several typical results and parametric studies 
are presented. Included here are cases run after the correlation 
studies (1) and following modifications to the VSAERO-TS program. 
The recent changes in the shedding model have significantly 
improved the treatment of the tapered and ogee planforms and 
additions to the wake influence coefficient evaluation for highly 
skewed panels have improved the treatment of the roll-up region 
near the tip vortex. A brief summary and conclusions are pre- 
sented in Section 8. 
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2.0 NOTATION 


B,C Velocity potential influence coefficients for constant 

source and doublet distributions, respectively, on 
quadrilateral panels, Eq. (7) 

c Airfoil chord (m) 

C p Pressure coefficient 

C L Lift Coefficient 

C M Pitching moment coefficient 

E Quantity in Eq. (7) 

f Frequency (Hertz) 

h Pitch axis unit vector 

i,j,k Orthogonal unit vector system defining the axes of the 
Cartesian coordinate system fixed relative to the 
blade. Figure 2 

£ Characteristic length, half mean chord (m) 

N Number of panels 

n Surface unit normal vector directed into the flow field 

R Position vector of a point relative to a point on the 

pitch axis 

S Blade surface 

s Surface distance (m) 

t Time (sec) 

T Normalized time, V t/ i 

00 

V Velocity (m/sec) 

v Perturbation velocity (m/sec) 

W Wake surface 

x,y,z Cartesian coordinates in the blade-fixed frame 
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a Onset flow incidence measured in the blade-fixed frame. 

Figure 2 

a Rate of change of a 

a Q Mean angle of attack 

Oscillation amplitude (mean to peak) 

A Finite increment 

n Spanwise location normalized by semispan 

y Doublet strength 

$ Total velocity potential 

<p Perturbation velocity potential 

a Source strength 

°R ,a M Source components due to unit pitch rotation and unit 

translation, respectively, Eq. (7) 

w Reduced Frequency 

V Gradient operator 

Subscripts 

J,K Values on panels, J,K, respectively 

U,L Upper, lower 

R,I Real, imaginary components normalized by a i 

W Wake 

00 Reference onset condition 
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3.0 MATHEMATICAL MODEL OVERVIEW 


3.1 FQLmulatign 

The general arrangement of a configuration is shown in Figure 
3 and is described in terms of a blade-fixed GLOBAL COORDINATE 
SYSTEM referred to as the G.C.S. The configuration is immersed 
in a uniform onset flow, V^, with velocity potential ^ . Poten- 
tial flow is assumed to exist both inside and outside the 
boundaries of the configuration. In the external flow field 
(i.e. f the one of interest) the total velocity potential, $ , is 
the sum of the onset flow potential, <p m , and the perturbation 
potential, 4> . Similarly, the total velocity potential inside the 
configuration is For the present description, the wake 
surfaces are assumed to have vanishing thickness with zero 
entrainment. 

After applying Green's Theorem to the inner and outer regions 
and combining the resulting expressions, the velocity potential 
at a point P on the inside surface can be written 



where r is the length of the vector from the surface element to 
the point, P, and S-P signifies that the point, P, is excluded 
from the surface integral. - $ L is the local jump in poten- 

tial across the wake surface, W. Equation (1) gives the total 
potential at the interior point, P, as the sum of perturbation 
potentials due to a normal doublet distribution of strength ($ - 
$i) on S and ($y - $ L ) on W, respectively, and a source distribu- 
tion of strength, n • ( “ V$) on S. The potential for the 

uniform onset flow, <P m , is also included. 

In principle, an infinite number of combinations of doublet 
and source distribution will give the same external flow field, 
but different internal flow fields. To render a unique combina- 
tion of singularities, either one of the singularity distribu- 
tions must be specified (e.g., a = 0 in the doublet-only formula- 
tion) or, as in the present case, the internal flow must be 
specified. There are several reasonable options for the internal 
potential flow. One of the most convenient options is to set 
= . This internal Dirichlet boundary condition gives zero 

perturbation to the onset flow inside the configuration. 
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Equation (1) becomes 


<J»n • V dS - 2 ir <pp 


($ 0 - $ L )n • V (i) dW 


^ n • V<f> (2) 


where <P, the perturbation potential in the flow field has been 
substituted for <J> - $ . 

co 

The first two terms in Eq. (2) give the perturbation poten- 
tial due to a distribution of normal doublets of strength, <f> , on 
the configuration surface, S. Similarly, the third term repre- 
sents a doublet distribution of strength, - $ L , on the wake 
and the fourth term represents a source distribution on strength 
a = n • V<|) on the configuration surface. 

For the problem of analyzing the flow about a given con- 
figuration geometry, the doublet distribution, , is the unknown 
while the source distribution is determined directly by the 
external Neumann boundary conditions; i.e., the source term in 
Eq. (2) can be evaluated directly from the condition of no flow 
penetration at the surface. The flow velocity relative to the 
blade-fixed frame is 



V = v + V ro - ah - R , (3) 

where the perturbation velocity, v = — V <f> . 


For zero penetration, V • n * 
n • V<p = n • V 

T oo 

and Eq. (2) becomes 



s 


0. Hence, 

- an • h -n R, 

+ II K - h) n * v (?) dw 

w 

- R ) dS (4) 
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This is the basic equation of the method. It is solved for 
the unknown surface perturbation potential, D, or surface doublet 
distribution, at a number of time steps as the blade proceeds 
through pitch oscillations. The wake surface is relocated at the 
end of each step to satisfy the zero load condition. The doublet 
distribution, 4>u “ • on the wa ke is related to conditions at 
the trailing edge in the steady case, while in the unsteady, 
time-stepping mode, the wake doublet values and wake location are 
known form solutions at earlier time steps. In the unsteady mode 
the Kutta condition 



3y 

9 


0 


(5) 


is satisfied at points along the wake separation line at each 
time step. 

At each time step the flow solution is determined with 
reference to the blade-fixed frame. The incompressible pressure 
coefficient is, therefore, given by 


C = (-V 2 - V 2 + 2 |£)/V 2 
P V s 9t/ 00 (6) 

where V g = ah ~ R - is the instantaneous velocity of a point 
on the surface relative to a stationary reference frame, and V is 
given by Eq. (3). 

Although the modulus of V x is constant in the present calcu- 
lations, its components in the blade-fixed frame are functions of 
time, with the instantaneous angle of attack 


a(t) = a + a. sin (pt) 

0 1 

Hence, the instantaneous rotation rate about the pitch axis is 

a (t) = ou p cos (pt) . 

In the case of harmonic (unsteady) flow, the motion and the 
velocities are assumed to undergo a periodic variation. For 

ipt iuiT 

example, the doublet strength, y = y Q = V e , where y is the 

amplitude, a reference length (usually semi-chord); u> = p^/V^ 
(reduced frequency); and T = tV^/ & (non-dimensional time). For 
harmonic flow, one solves for the amplitude of doublet strength 
(complex), taking into account its variation in the wake con- 
sistent with the unsteady Kutta condition. 
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4.0 NUMERICAL PROCEDURE 


The numerical procedure to obtain the unsteady pressure 
distribution, forces and moments is shown in Figure 4. The 
surface of the blade is represented by planar quadrilateral 
panels over each of which the doublet and source distributions 
are assumed constant. The input data is assembled in three main 
parts: basic data, surface patch geometry description and basic 
wake geometry description. Also, automatic paneling routines are 
installed to simplify user input. Complete details of the input 
parameters, output description, a detailed diagram showing the 
connections between subroutines in the code and user notes for 
the geometry description are presented in Reference 8. 


Equation (4) is satisfied simultaneously at a point at the 
center of each panel. If there are N panels representing the 
blade surface, Eq. (4) becomes: 


N 


Z K 


KV 



2it ^j + e j = 0; 


J=1,N 


(7) 


where 

4> k /4tt . ) K 


is the unknown doublet value on panel K. 


n T7 

w 

Z 



+ ota 


R. 



(Note: U K = 


where N w , the number of panels in the wake, varies with time and 

a and V take their instantaneous values at each time step. 

00 


a 




B 


JK 


/ 4 TT 


is the source distribution due to blade rotation about the pitch 
axis, and 



n 


K 



/ 4 TT 


are the components of a three-part source distribution due to 
the relative translation of the blade and the onset flow. (Note: 
in the present symmetrical case the y-component is zero.) 
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Figure 4. Flow Diagram for the Method. 
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The quantities, Bjk and C jk* are the velocity potential 
influence coefficients for the constant source and doublet dis- 
tributions, respectively, on panel K acting on the control point 
on panel J. These include contributions from the image panel in 
this case. Expressions for these influence coefficients have 
been given by Morino (9) based on hyperbolic paraboloidal panels. 
Slightly different expressions are installed in the present code 
based on planar panels. 

When the influence coefficient matrix has been assembled, 
Eq. (7) is solved by a direct method for N ^ 320 and by an 
iterative method for N > 320. 

The surface pressure distribution is calculated using Eq. 
(6). The surface gradient of y is evaluated on each panel by 
differentiating a two-way parabolic fit through the doublet 
values on the panel and its four immediate neighbors. At the 
separation lines a simple differencing is applied for the 
gradients approaching the separation line. 

The gradient of $ with respect to time is evaluated by 
central differencing over two time steps; i.e., 

3d) / t t-2At\ 

If = U - * j/2At . 


The real and imaginary pressures are obtained by Fourier 
analysis for the first harmonic based on solutions over a com- 
plete cycle. The calculations start with incidence a Q and a 
regular (i.e., steady) wake. Two iterations are performed to 
render the wake force free. Based on several test cases, two 
iterations proved to be adequate. An oscillatory doublet com- 
ponent based on a linearized solution is then superimposed along 
each wake line before starting the time-step mode. Time-step 
calculations proceed over a half cycle before applying the 
Fourier analysis. This was done to eliminate the transient 
effects of the initial calculation. 


At each time step a new panel is formed at the head of each 
column of wake panels and all the existing wake panel corner 
points are convected downstream at the computed local velocity. 
The doublet value on each new panel is based on the conditions at 
the separation line and satisfies the Kutta condition, Eq. (5). 
It is assumed that the shedding occurs at constant vorticity over 
the time interval. At. In this way the doublet strength, y w t+At , 
on the new wake panel is related to the strength, y , of the 
previous wake panel at the trailing edge by 


where 



t+At 



y 


t 


w 


is the resultant doublet value at the separation line. 
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Preliminary tests of the time-stepping calculation showed 
some sensitivity to time-step size relative to the local panel 
size at the trailing edge on each strip. In order to alleviate 
this a scheme was developed in the program whereby the x-wise 
location of the wake panels were fixed for all time by the user. 
Thus, at the outset of a calculation the user specifies a set of 
x-stations as locations of vertical cross-flow planes which will 
contain wake panel corner points. Figure 5. In principle, the 
planes can be located with a reasonable relationship to wake 
shedding panels at the blade trailing edge. The plane spacing 
can be made dense near the blade and can be opened out with 
increasing distance downstream. At each time step points on the 
present wake are first transported a small distance according to 
the local velocity and a time-step size; the NUWAKE routine then 
interpolates (using a biquadratic scheme) along each line and 
reforms the wake panel corner points — and the wake panel doublet 
values — in the user-specified wake-grid-plane system. Figure 6. 
This scheme alleviates the sensitivity to time-step size and also 
leads to a more effective use of a limited number of wake panels. 
Initially (1) the NUWAKE model had some difficulties in treating 
the tapered and ogee planforms, especially at the smaller reduced 
frequencies. The scheme has now been modified to give an im- 
proved trailing-edge doublet value (using an extrapolation from 
values at previous times) for interpolation in the immediate zone 
near the blade trailing edge. This has given a significant 
improvement in the behavior of the method relative to that in 
Reference 1, and good solutions are now obtained for the tapered 
tip and ogee planforms (Section 7). 

In the preliminary version of the time-stepping program the 
wake panels had a piecewise constant doublet distribution. Later 
this was changed to a linear doublet distribution (i.e., constant 
vorticity) in the streamwise direction to provide a smoother 
representation of the vorticity. The influence coefficient for 
each panel is formulated for a planar surface and so a skewed 
panel is represented by a projected flat quadrilateral lying in 
the mean plane. Panels in the extreme roll-up region of the tip 
vortex are highly skewed and consequently the flat-panel influ- 
ence coefficient representation leaves a large number of "holes" 
in the apparent vortex sheet surface. These holes cause a loss 
in the local induced effects on the blade tip surface and this 
leads to errors in the computed pressure distribution near to the 
tip vortex. The program has therefore been modified to treat 
each highly skewed panel as a pair of triangles. Figure 7. One 
of the diagonals of the quadrilateral is taken as the common side 
for the two triangles. In this way the influence coefficient 
evaluation now covers a continuous surface since the edges of 
skewed panels match up with their neighbors. The panels are 
still regarded as quadrilaterals as far as geometry, etc. is 
concerned; the triangle treatment is only applied in the near- 
field influence routine. 
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Figure 7. Treatment of Highly Skewed Wake Panels. 
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5.0 1 VSAERO-TS ' : VISCOUS/POTENTIAL FLOW ITERATIVE SCHEME 


The time-stepping version of the code includes a viscous/po- 
tential flow iteration loop which is used to establish the extent 
of tip-edge separation. The calculation is performed in the 
steady mode at the mean angle of attack, the results of which can 
be used to set the extent of tip-edge vortex separation for the 
unsteady calculation. 

The rectangular tip was used to test the procedure for 
calculating the location of the tip-edge separation and boundary 
layer calculations. Presented here is a sample case. 

Following a steady calculation at the mean angle of attack, 
a family of surface streamlines were calculated. Figure 8. 
Boundary layer layer calculations were then performed along each 
streamline. The boundary layer code is an integral method based 
on locally axisymmetric conditions with the surface curvature and 
streamline convergence/divergence effects included. First, cal- 
culations were performed in the steady mode to establish the 
extent and location of separation lines as a function of angle of 
attack. Angles of 4°, 8°, 12° and 16° were analyzed. In each 
case a family of streamlines was calculated. Figure 9 shows the 
extent of separation as a function of angle of attack and two 
Reynolds numbers. The lines represent the conditions for the 
first iteration; i.e., before modeling the edge vortex. The low 
Reynolds number corresponds with the unsteady experiment while 
the higher Reynolds number is more representative of full-scale 
conditions. 

Figure 10 shows the plan view of the calculated separation 
boundaries for a range of angles of attack. In addition to the 
tip-edge separation, the separation boundaries over the inner 
part of the wing are also included. Figure 10(a) shows the 
boundaries for the model Reynolds number of 1.28 x 10 6 and Figure 
10(b) is for a full-scale Reynolds number of 6.0 x 10. At the 
lower Reynolds number, some laminar separation with no reat- 
tachment is indicated at 18° in the inboard region. At the 
higher Reynolds number, the extent of separation is considerably 
lower and would probably have some impact on dynamic stall quali- 
ties relative to the model data. For the purposes of the model 
scale Reynolds number, the inboard separation appears to be small 
enough to be neglected for angles of attack up to 120. 
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Figure 8. Calculated Streamlines on Rectangular Tip. 
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6.0 1 VSAERO-H ' : DISCUSSION ON WAKE INTEGRAL AND CONVERGENCE 

CHARACTERISTICS 


Program VSAERO-H uses harmonic wake and small amplitude 
assumptions. The governing flow equation is 


N 


s 


2> 

Tf=1 


K 


K^J 


C ) 
JK 



K=1 


(y c ) 

WK JK 


+ y: { °K b jk ) 

K=1 


0 ; 


J=l, N s 


( 8 ) 


where N s and N w are the number of surface and wake panels. 

In the case of harmonic flow f the application of the un- 
steady Kutta condition enables one to express the doublet 
strength in a wake strip in terms of its strength at the begin- 
ning of the strip. The wake contribution in Eq. (8) is replaced 


yi y WK C JK' 
K=1 


where u Wg is the number of spanwise (wake) strips and U WK is the 
doublet strength at the beginning of the k th wake strip, and 


C 


JK 



-iw (x-x )n 
e s 


Wake Strip K 


V 


1 


r 


J 


dS. 


The convergence characteristics of the wake integral were 
thoroughly investigated by adopting the following mathematical 
model. The wake influence coefficients (40 are computed using an 
analytic expression for spanwise integration and a numerical 
scheme for chordwise integration from the trailing edge to in- 
finity. However, due to the highly convergent nature of the 4 1 
integrand, which is inversely proportional to r 3 , it is only 
necessary to integrate a few chord lengths behind the trailing 
edge. A series of convergence studies are conducted to evaluate 
the computational efficiency of the wake influence coefficients. 
Two singularity models, a constant doublet distribution and a 
linear doublet distribution over each wake panel were evaluated. 
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Figure 11 shows the model used for the numerical scheme in 
the present study. The doublet distribution in the wake is 
assumed to be 


y (1 - y/b) 2 e~ ia)X 

where y Q is the doublet strength at the mid-span trailing edge of 
the wing, b is the semi-span, and w is the reduced frequency. 
The wake integrals are evaluated by summing the contributions of 
several spanwise strips covering the entire wing span. For 
chordwise integration, two doublet models, a constant distribu- 
tion and a linear distribution are used. For the constant doub- 
let distribution model, the doublet strength at any spanwise 
strip. 


y 


t 


-iwx 

e 


is assumed to be a constant and is equal to the strength at the 
center of the chosen wake strip. For example, the doublet 
strength at the n th wake strip is assumed to be 

-iu3X n 


where x n is the center of the n wake strip, and the contribu- 
tion is computed numerically over the strip width, Ax n . For the 
linear doublet distribution model, a wake strip influence is 
computed by summing the contribution of the two terms, the doub- 
let strength at the leading edge of the strip. 


-ito (x - x / 2 ) 

y t e 

and the incremental doublet strength, 

-i u x_ -io) 


y t 


n 


- e 


io) (x n ~Ax n /2) 


Figure 12 shows the convergence characteristics of the wake 
influence coefficients as the panel size in the wake varied from 
1/4 to 4 widths of the trailing-edge panel. As can be seen from 
Figure 11, the linear doublet distribution model improves the 
convergence characteristics and hence this model is used in the 
program. 

Additionally the numerical results for <p were compared with 
those of analytical results presented in Reference 10. The wake 
model used for this computation is presented in Figure 13 and the 
comparison of the results for both wake models, constant and 
linear doublet, are presented in Table I. As can be seen, the 
comparison is excellent, validating the numerical scheme adopted 
in the current program. 
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Figure 11. Model used for Convergence Study. 
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Figure 13. Model used for Wake Integral Computations. 
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CO 

X 

ANALYTICAL 

(Send) 

COMPUTED 

Ax = .025 

Ax = 

.1 

Ax = .4 

.1 

-1.5 

1.233 - 

.2861 

1.236 - . 286i 

1.235 - 

. 286i 

1.214 - . 288i 

.1 

-1.0 

2.810 - 

. 472i 

2.827 - .473i 

2.822 - 

.4731 

2.727 - . 480i 

.1 

-.5 

10.300 - 

1.012i 10.502 - 1.022i 

10.428 - 

1.026i 

9.454 - 1.057i 

.5 

-1.5 

.700 - 

. 551i 

.700 - . 554i 

.699 - 

.5551 

.678 - .5661 

• 5 

i 

• 

O 

1.915 - 

1. 140i 

1.927 - 1. 114i 

1.921 - 

1.115i 

1.825 - 1. 179i 

.5 

-.5 

8.510 - 

3 . 221i 

8.701 - 3 . 257i 

8.627 - 

3.27H 

7.644 - 3.426i 


LINEAR DOUBLET MODEL; $ r X 10 3 t Y = ♦ 1 , Z = .1 . 


(0 

X 

ANALYTICAL 

(Send) 

COMPUTED 

Ax = 

.025 

Ax 

- 

.i 

Ax = .4 

.1 

-1.5 

1.233 - . 286i 

1.237 - 

. 286i 

1.237 

- 

. 286i 

1.236 - . 289i | 

.1 

-1.0 

2.810 - .472i 

2.827 - 

.473i 

2.827 

- 

. 474i 

2.827 - .484i 

.1 

-.5 

10.300 - 1.012i 10.507 - 

l.OlOi 

10.507 

- 

1. 028i 

10.505 - 1.073i 

.1 

-1.5 

.700 - . 551i 

.700 - 

.5541 

.700 

- 

.555i 

.697 - . 571i 

.1 

-1.0 

1.915 - 1. 140i 

1.923 - 

1 . 144i 

1.927 

- 

1.1491 

1.918 - 1.1941 

.1 

-.5 

8.510 - 3.221i 

8.706 - 

3 . 257i 

8.703 

— 

3.2821 

8.667 - 3.498i 


Table 1. Comparison of Computed and Analytical Unsteady 
Velocity Potential ($p) . 
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7.0 RESULTS AND DISCUSSION ON CONVERGENCE CHARACTERISTICS 


In this section a limited number of results are presented 
for the four tip planforms to demonstrate the convergence charac- 
teristics of VSAERO-TS as a function of reduced frequency and the 
number of time steps per cycle. A complete set of the results 
may be found in Reference 1 in comparison with experimental data. 


7.1 Rectangular Tip — Parametric Studies 

An extensive set of calculations have been performed using 
VSAERO-TS. Basically, calculations were performed using 60 to 
120 time steps over 3 half cycles. 

(1) E£f,ect..Qf Panel Density 

First, the effect of panel density was examined: 

(a) 30 (chordwise) x 4 (spanwise) plus 3 (vertical) x 
15 (chordwise) on the tip 

(b) 40 x 4 + 3 x 20 

(c) 30 x 8 + 3 x 15 


The calculations were performed for to = 0.3 (based on half- 
chord), otj_ = 1° and a 0 = 12°. The wake length was 10 chords with 
25 wake panels down each strip and N = 120 (time steps). The 
effect of chordwise panel density is given in the table below 
which compares the unsteady lift and moment coefficients at two 
spanwise stations. 


Spanwise Panels 

Station 

y/s 30 x 4 + 3 x 15 40 x 4 + 3 x 20 


0.25 



4.42 + 1. 17i 
0.573 - 0 . 517i 


4.42 + 1.171 
0.563 - 0 . 513i 


0.95 



2.52 + 0 . 86i 
0.098 - 0 . 395i 


2.52 + 0 . 86i 
0.079 - 0395i 


Table 2. Effect of Chordwise Panel Density on Spanwise Lift and 
Moment Distributions. 
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The results are therefore essentially converged with 30 
chordwise panels. In fact, they are converged even for as few as 
20 chordwise panels. The chordwise distribution of panels is the 
default form, which increases panel density towards the leading 
edge. 


The effect of spanwise panel density is shown in Figure 14 
for the unsteady lift distribution. Four spanwise intervals 
clearly give a very good account; however, for details of the 
behavior near the tip the eight spanwise intervals case probably 
represents a minimum density. Both distributions used a "cosine" 
spacing with density increasing towards the tip. 


(2) Effect o f Wake Length 

The effect of wake length (i.e., the length of wake repre- 
sented by the wake doublet panels) was examined. Presented 
below are the integrated unsteady lift and moment characteristics 
at the inboard and outboard stations, for the case of o =12°, 
= 1°, and u> = 0.30. 


Spanwise Wake Length 

Station 

y/s 10 Chords 20 Chords 


0.25 


0.95 


C L 

4.48 

+ 

1.24i 

4.47 

+ 

1. 25i 

C M 

0.038 

- 

0 . 441i 

0.042 

- 

0 . 439i 

C L 

2.27 

+ 

0 . 88i 

2.27 

+ 

0. 88i 

% 

-0.011 

- 

0.280i 

-0.011 


0 . 280i 


Table 3. Effect of Wake Length on Spanwise Lift and Moment Dis- 
tributions. 


The results are essentially converged with a wake length of 
10 chords. Although not presented in this document, the wake 
length of 10 chords proved to be adequate even for the case of 
low reduced frequency of 0.1. 
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Figure 14. Calculated Unsteady Lift Distribution--Ef feet 
of Spanwise Panel Density (a = 12°, a. = 1°, 
aj = 0.30) . ° 1 



(3) Effect of Number, of Time Steps 


For the case of a 0 = 12°, “i = 1°, id = 0.30, and wake length 
of 10 chords, the number of time steps was varied between 40 and 
120. Table 4 compares the integrated unsteady lift and moment 
coefficients at several spanwise stations. 

This example demonstrates that N = 60 is adequate to obtain 
convergent results. However, for the low reduced frequencies, N 
= 120 is required for the proper convergence. 

Additionally, some typical results for the rectangular tip 
are presented in Figures 15 and 16 for reduced frequencies of 0.1 
and 0.3, respectively. The chordwise pressure distributions at 
two spanwise stations (y/s and 0.25 and 0.95) are compared with 
those of the DFVLR test results. Figures 15(a) and (b) show the 
chordwise distributions at the two spanwise stations for w = 0.1; 
the computed results for these are obtained using the standard 
VSAERO-TS quadrilateral panels. For the present case of high 
angle of attack, this approach is inadequate to handle the roll- 
up of the tip vortex due to the problem of panel skew. Hence, 
the modified VSAERO-TS, in which the quadrilateral panels in the 
wake are temporarily replaced by pairs of triangular panels, is 
used and the results are shown in Figures 15(c). As can be seen 
the comparison between the computed and test results has substan- 
tially improved. Figures 15(d) and (e) show the instantaneous 
tip vortex wake geometry for this particular case. Figures 16(a) 
and (b) show the chordwise pressure distributions for w = 0.3 at 
the root and tip sections. In general, the comparison is good 
for both real and imaginary parts. 


7.2 Swept Tip 

This particular case is used to investigate the effect of 
wake transport length for each time step. In order to maintain 
the same level of wake detail, the number of time steps for each 
cycle in VSAERO-TS must be inversely proportional to the reduced 
frequency. The swept tip case is considered and the number of 
time steps are varied between 120 and 40 as reduced frequency is 
varied between 0.1 and 0.3. 

Figures 17(a), (b) and (c) show the comparison of the 
chordwise pressure distribution between VSAERO-TS and the DFVLR 
test results at a 60% spanwise location for ui - 0.1, 0.2 and 0.3, 
respectively. As can be seen, the comparison is uniformly good 
in all these cases. 
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Table 4. Effect of Number of Time Steps on Spanwise Lift and Moment Distributions. 
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Figure 15 (b) . Comparison of Chordwise Pressure Distribution at 
y/s = 0.95 between Computed (VSAERO-TS) and DFVLR 
Test, Rectangular Tip (a = 12° , a. = 1.066° , w = 
0 . 1 ). ° 1 
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Figure 15(c). Comparison of Chordwise Pressure Distribution at 

y/s = 0.95 between Computed (VSAERO-TS, Triangular 
Panels) and DFVLR Test, Rectangular Tip (a = 12° , 
a. = 1.066° , a) = 0.1) . 



Figure 15(d). Instantaneous Wake Geometry — Plan View. 
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0*0 0.2 0.4 0.6 0.8 1.0 


x/c 



0.0 0.2 0.4 0.6 0.8 1.0 

x/c 


(b) At y/ s = 0.60 (a Q — 4° , — 0.710° , to — 0.20, N 80) . 

Figure 17. Continued. 
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x/c 



0*0 0.2 0.4 0.6 0.8 1.0 

x/c 

(c) At y/s = 0.60 (a Q = 4° , a ± = 0.710° , to = 0.30, N = 40) . 

Figure 17. Concluded. 
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7.3 Oge e . lip 


One case of comparison of the ogee tip chordwise pressure 
distribution between VSAERO-TS, VSAERO-H and the DFVLR test re- 
sults is presented in Figures 18(a) and (b) for w = 0.3 at an 85% 
spanwise location. For this particular case the harmonic wake 
results appear to correlate better with the test results. How- 
ever, boundary layer displacement effect is not modelled in the 
analysis at this time and would be expected to reduce the suction 
levels in both computed pressure distributions. 


7.4 Spanwise Distribu tion of Unsteady. Lift — Four Blade Tips 

Figures 19(a) and (b) compare the spanwise distributions of 
real and imaginary lift coefficients for the four planforms for to 
= 0.1 and 0.3, respectively. It is emphasized that these coeffi- 
cients are based on local chord. The swept and rectangular tips 
have almost the same distributions. The slightly larger imagin- 
ary part on the swept tip is consistent with the heaving effect 
caused by the rearward location of the local quarter-chord rela- 
tive to the pitch axis. 

A similar heaving effect is present for the tapered tip 
except in this case the decreasing chord lowers the effective 
value of the reduced frequency. Because of this, the imaginary 
term has taken negative values near the tip. The real part of 
the lift component on the tapered tip take increased values 
towards the tip because of the reducing chord; this is consistent 
with steady solutions. This effect is even more pronounced on 
the ogee tip; here, the increased real lift values would probably 
lead to separations in the outboard region. The large negative 
movement in the imaginary component near the ogee tip is consis- 
tent with (a) the small local chord lowering the effective value 
of the reduced frequency; and (b) the heaving effect caused by 
the local quarter-chord points located ahead of the pitch axis. 
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(b) Imaginary Part. 
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Figure 18. Concluded. 
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Figure 19 . Calculated Spanwise Distributions of Unsteady 
Lift for the Four Planforms (VSAERO-TS) . 
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(b) a Q = 12° , a. = 1.066° , a> = 0.3. 
| Figure 19. Concluded. 
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8.0 SUMMARY AND CONCLUSIONS 


In this report the theoretical development of two computer 
programs, VSAERO-TS and VSAERO-H, for calculating the unsteady 
aerodynamic characteristics of arbitrarily-shaped wings oscil- 
lating in pitch, is presented. The effect of several parame- 
ters — wake length, the number of time steps, and panel density — 
on the chordwise pressure distribution in VSAERO-TS is presented. 
Also, the convergence characteristics of both programs are dis- 
cussed. 

The treatment of the unsteady Kutta condition has been 
improved by a refined wake shedding model. Whereas the old model 
would not give consistent solutions for the ogee and tapered tip 
planforms, the new treatment in VSAERO-TS gives good solutions 
for these planforms. 

In the computer programs the influence coefficient for each 
panel is formulated for a planar surface and so a skewed panel is 
represented by a projected flat quadrilateral lying in the mean 
plane. Panels in the extreme roll-up region of the tip vortex 
are highly skewed; for the proper representation of these cases, 
the program has been modified to treat each highly skewed panel 
as a pair of triangles. 

The two programs are validated by comparing the chordwise 
pressure distribution of several blade tip planforms with experi- 
mental data. The comparison, for the most part, is good. The 
triangular panel representation improved the chordwise pressure 
distribution near the tip region for higher mean angle of attack. 
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